﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace _708DllTest
{
    class GeLinNiZhi
    {

        /// <summary>
        /// 半长轴
        /// </summary>
        public double BanChangZhou;
        /// <summary>
        /// 偏心率
        /// </summary>
        public double PianXinLv;
        /// <summary>
        /// 轨道倾角
        /// </summary>
        public double GuiDaoQingJiao;
        /// <summary>
        /// 升交点赤径
        /// </summary>
        public double ShengJiaoDianChiJing;
        /// <summary>
        /// 近地点角
        /// </summary>
        public double JinDiDianJiao;
        /// <summary>
        /// 周期
        /// </summary>
        public double ZhouQi;
        /// <summary>
        /// 开始纪元时刻的平近点角
        /// </summary>
        public double PingJinDianJiao;


        /// <summary>
        /// 开始纪元 格式："2018 04 18 10:15:00"或"2018-04-18 10:15:00"
        /// </summary>
        public string KaiShiJiYuan;
        /// <summary>
        /// 开始纪元的时间戳
        /// </summary>
        private double StartTime;


        private string Msg;




        /// <summary>
        /// 角度转弧度值
        /// </summary>
        /// <param name="angel">角度</param>
        /// <returns>弧度</returns>
        private double AngleToHudu(double angle)
        {
            return Math.PI * angle / 180.0;

        }

        /// <summary>
        /// 计算卫星周期
        /// </summary>
        /// <param name="a">半长轴</param>
        /// <returns></returns>
        public double Com_ZhouQi(double a)
        {
            return (double)2 * Math.PI * Math.Sqrt(Math.Pow(a / 1000, (double)3) / 398600.44);
        }





        // Update is called once per frame
        public string Update()
        {



            //计算的时间
            string time = "2018-01-01 12:00:00";// 例如： time  = 2018-01-01 12:00:00



            //NewComplete(BanChangZhou, PianXinLv, GuiDaoQingJiao, ShengJiaoDianChiJing, JinDiDianJiao, time);
            NewComplete(42164142, 0, 55, 338, 0, time);

            return Msg;

        }

        public List<double> NewComplete(double a, double e, double i, double sjdcj, double w, string nowtime)
        {
            List<double> JWD = new List<double>();

            //获取当前时间的年月日时分秒

            DateTime dt = Convert.ToDateTime(nowtime);
            int Year = dt.Year;
            int Mon = dt.Month;
            double Day = dt.Day;
            double st = dt.Hour * 3600 + dt.Minute * 60 + dt.Second;


            //*************************** 第一步   求当前时间的平近点角（这个部分不确定是否要转换到0-360之间，测试发现转与不转 最后计算的f值都是一样的）  ***********************************
            double NowTime = (Convert.ToDateTime(nowtime).ToUniversalTime().Ticks - 621355968000000000) / 10000000;


            ZhouQi = Com_ZhouQi(a);


            //StartTime = (Convert.ToDateTime("2000-01-01 00:00:00").ToUniversalTime().Ticks();
            //当前时间所转过的平近点角度
            double angle = PingJinDianJiao + ((NowTime - StartTime) / ZhouQi) * 360;
            //将角度转换到0-360之间
            double M = (angle / 360 - Math.Floor(angle / 360f)) * 360;

            //double M = PingJinDianJiao + ((NowTime - StartTime) / ZhouQi) * 360;

            //*************************** 第二步   平近点角—>偏近点角->真近点角 ***********************************
            double E = M2E(e, M);

            double f = E2f(e, E);

            //*************************** 第三步   计算赤经和赤纬  ***********************************

            double u = w + 211.500; //卫星离升交点的角距

            double cj = sjdcj + Math.Atan(Math.Tan(AngleToHudu(u)) * Math.Cos(AngleToHudu(i))) / Math.PI * 180.0; //计算赤经

            double cw = Math.Asin(Math.Sin(AngleToHudu(u)) * Math.Sin(AngleToHudu(i))) / Math.PI * 180.0;  //计算赤纬

            // Debug.Log("赤经：" + cj);
            //*************************** 第四步   计算格林威治的恒星时角得出地心经纬度 ***********************************

            double jd = cj - calGAST(Year, Mon, Day, st);
            double wd = cw;




            //Debug.Log("半长轴：" + a + "  偏心率：" + e + " 轨道倾角：" + i + "  升交点赤经：" + j + "  近地点角：" + w +"平近点角："+PingJinDianJiao);

            Msg += "经度：" + jd + "  纬度：" + wd + "    周期：" + ZhouQi;
            //Debug.Log("卫星：" + gameObject.name + "经度：" + jd + "  纬度：" + wd + "    周期：" + ZhouQi);
            Msg += "   时间：" + nowtime;
            //Debug.Log("格林尼治恒星时角：" + calGAST(Year, Mon, Day, st));
            //Debug.Log("时间string:" + nowtime + " " + Year + " 年 " + Mon + " 月  " + Day + "日  st:" + st);
            // Debug.Log("angel:" + +angle/360 + "  平近点角:" + M);
            // Debug.Log(" 平近点角:" + M);
            // Debug.Log("*********************************************************************************************");

            JWD.Add(jd);
            JWD.Add(wd);
            return JWD;

        }



        #region 平近点角、偏近点角、真近点角装换

        /// <summary>
        /// 平近点角->偏近点角
        /// </summary>
        /// <param name="e">偏心率</param>
        /// <param name="meanAnom">平近点角（度数）</param>
        /// <returns>偏近点角（度数）</returns>
        private double M2E(double e, double meanAnom)
        {
            double eccAnomX = 0.0;
            meanAnom = meanAnom * Math.PI / 180.0;//度数转化弧度
            double e3 = e * e * e;
            double e2 = e * e;
            eccAnomX = meanAnom + (e - (e3 / 8.0)) * Math.Sin(meanAnom) + 1.0 / 2.0 * e2 * Math.Sin(2.0 * meanAnom) + 3.0 / 8.0 * e3 * Math.Sin(3.0 * meanAnom);
            eccAnomX = eccAnomX * 180.0 / Math.PI;//弧度转化为度数
            return eccAnomX;
        }
        /// <summary>
        /// 偏近点角->真近点角
        /// </summary>
        /// <param name="e">偏心率</param>
        /// <param name="E">偏近点角（度数）</param>
        /// <returns>真近点角（度数）</returns>
        private double E2f(double e, double E)
        {
            double nuX = 0;
            E = E * Math.PI / 180.0;//度数转化成弧度
            nuX = 2.0 * Math.Atan(Math.Sqrt((1.0 + e) / (1.0 - e)) * Math.Tan(0.5 * E));
            nuX = nuX * 180.0 / Math.PI;//弧度转化成度数
            return nuX;
        }

        #endregion



        #region 求格林威治的恒星时角

        const long IGREG = (15 + 31 * (10 + 12 * 1582));

        /// <summary>
        /// 计算儒略日
        /// </summary>
        /// <param name="Y">年</param>
        /// <param name="M">月</param>
        /// <param name="D">日</param>
        /// <returns>给定时刻的儒略日</returns>
        private int GetJD(int Y, int M, double D)
        {
            int jd;

            int A, jy = Y, jm;

            if (jy == 0)
            {
                //Debug.Log("错误！年份为零！");
                return 0;
            }
            if (jy < 0)
            {
                ++jy;
            }

            if (M > 2)
            {
                jm = M + 1;
            }
            else
            {
                --jy;
                jm = M + 13;
            }

            jd = (int)(Math.Floor(365.25 * (jy + 4716)) + Math.Floor(30.6001 * jm) + D - 1524.5);
            if (D + 31 * (M + 12 * Y) >= IGREG)//判断是否为格里高利历日IGREG  
            {
                A = (int)(0.01 * jy);
                jd += 2 - A + (int)(0.25 * A);//加百年闰  
            }
            return jd;
        }

        /// <summary>
        /// 计算任意时刻的格林尼治平恒星时
        /// </summary>
        /// <param name="jd">儒略日</param>
        /// <param name="st">当前时间秒（时*3600+ 分*60 + 秒）</param>
        /// <returns></returns>
        private double calGMST(double jd, double st)
        {
            double D = jd - 2451545.0;//Compute the number of days since Jully2000  
            double T = D / 36525.0;
            double GMST = 6.697374558 + 2400.051336 * T + 0.000025862 * T * T + st / 3600 * 1.00273790935;
            GMST = GMST % 24 * 15;
            return GMST;
        }
        /// <summary>
        /// 计算任意时刻的格林尼治恒星时角
        /// </summary>
        /// <param name="yr">年</param>
        /// <param name="mo">月</param>
        /// <param name="day">日</param>
        /// <param name="st">当前时间秒（时*3600+ 分*60 + 秒）</param>
        /// <returns>格林尼治恒星时角</returns>
        private double calGAST(int yr, int mo, double day, double st)
        {
            double fjd = GetJD(yr, mo, day);
            double TJD = fjd - 2451545.0;//2451545为2000 1月1.5的儒略日  
            double T0 = TJD / 36525.0;
            double THETAm = calGMST(fjd, st);
            double EPSILONm = 23.4392911111111111 - 0.0130041666667 * T0 - 1.638e-07 * T0 * T0 + 5.0361e-07 * T0 * T0 * T0;

            double L = 280.4665 + 36000.7698 * T0;
            double dL = 218.3165 + 481267.8813 * T0;
            double OMEGA = 125.04452 - 1934.136261 * T0;
            double dPSI = -17.20 * Math.Sin(OMEGA) - 1.32 * Math.Sin(2 * L) - 0.23 * Math.Sin(2 * dL) + 0.21 * Math.Sin(2 * OMEGA);
            double dEPSILON = 9.20 * Math.Cos(OMEGA) + 0.57 * Math.Cos(2 * L) + 0.10 * Math.Cos(2 * dL) - 0.09 * Math.Cos(2 * OMEGA);

            //Convert the units from arc-seconds to degrees  
            dPSI /= 3600;
            dEPSILON /= 3600;

            double GAST = (THETAm + dPSI * Math.Cos(EPSILONm + dEPSILON)) % 360;
            return GAST;

        }




        #endregion
    }
}
